home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / randist / bigauss.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.0 KB  |  70 lines
  1. /* randist/bigauss.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_rng.h>
  24. #include <gsl/gsl_randist.h>
  25.  
  26. /* The Bivariate Gaussian probability distribution is 
  27.  
  28.    p(x,y) dxdy = (1/(2 pi sigma_x sigma_y sqrt(r))) 
  29.                     exp(-(x^2 + y^2 - 2 r x y)/(2c)) dxdy     
  30.  
  31. */
  32.  
  33. void
  34. gsl_ran_bivariate_gaussian (const gsl_rng * r, 
  35.                 double sigma_x, double sigma_y, double rho,
  36.                 double *x, double *y)
  37. {
  38.   double u, v, r2, scale;
  39.  
  40.   do
  41.     {
  42.       /* choose x,y in uniform square (-1,-1) to (+1,+1) */
  43.  
  44.       u = -1 + 2 * gsl_rng_uniform (r);
  45.       v = -1 + 2 * gsl_rng_uniform (r);
  46.  
  47.       /* see if it is in the unit circle */
  48.       r2 = u * u + v * v;
  49.     }
  50.   while (r2 > 1.0 || r2 == 0);
  51.  
  52.   scale = sqrt (-2.0 * log (r2) / r2);
  53.  
  54.   *x = sigma_x * u * scale;
  55.   *y = sigma_y * (rho * u + sqrt(1 - rho*rho) * v) * scale;
  56. }
  57.  
  58. double
  59. gsl_ran_bivariate_gaussian_pdf (const double x, const double y, 
  60.                 const double sigma_x, const double sigma_y,
  61.                 const double rho)
  62. {
  63.   double u = x / sigma_x ;
  64.   double v = y / sigma_y ;
  65.   double c = 1 - rho*rho ;
  66.   double p = (1 / (2 * M_PI * sigma_x * sigma_y * sqrt(c))) 
  67.     * exp (-(u * u - 2 * rho * u * v + v * v) / (2 * c));
  68.   return p;
  69. }
  70.